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ABSTRACT 


The  use  of.  the  Seqxisntial  Regression  Algorithm  (SER)  to 
coherently  remove  background  noise  from  an  ELF  sensor  is  presented. 
The  SER  algorithm,  is  described  for  a  multi-channel  application  in 
order  to  cancel  coherent  portions  of  reference  sensors  from  a 
primary  sensor.  The  algorithm  adaptively  accouiits  for  differences 
betv/een  two  parallel  array  platforms  for  the  purpose  of  coherent 
subtraction.  A  section  on  likelihood  detector  schemes  is  also 
presented.  This  v.’ork  is  in  support  of  a  submerged  ELF  sensor  array 
project  run  by  the  Jchr'.s  Hopkins  University  Applied  Physics  Lab. 
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I .  INTRODUCTION 


A.  BACKGROUND 

Extremely  low  frequency  (ELF)  geoelectric  noise  has  a 
highly  impulsive  nature  v/hich  maizes  it  very  difficult  t: 
remove  by  filtering.  The  main  contributors  to  che 
impulsiveness  of  the  ELF  spectrum  are  lightning  discharges. 
The  lov;est  frequencies  of  the  EL"  range  (1  to  60Hzi  have  such 
a  lev;  atte.nuation  rate  as  they  propagate,  that  iight:na 
discharges  ail  over  the  world  must  be  considered  v/hen  studying 
the  noise  spectrum  at  a  single  location.  ELF  signals  have 
been  stulied  for  such  applications  as  power  transfer, 
comrr.unications,  geophysical  surveys,  and  so  forth.  The  amount 
of  information  on  the  subject  is  immense  and  increasing.  How 
to  account  for  the  impulsive  background  electrom^agnetic  fields 
inevitably  present  in  any  system,  operating  in  this  frequency 
range  remains  a  serious  problem,  that  any  system  designer  must 
answer  . 


The  background  noise  of  the  lower  ELF  frequency  range  can 
be  thought  of  as  consisting  of  both  local  and  gicbal 
phenom.ena .  The  global  component  is  comprised  of  signals  that 
are  relatively  coherent  over  extended  distances.  A  detailed 
recording  of  data  in  this  frequency  range 

...is  characterized  by  bursts  of  a  few  cycles  duration 
which  miay  be  ten  or  more  times  the  mean  amplitude  between 


1 


bursts,  This  bursty  character  is  apparently  coherent  over 
most,  if  not  all,  of  the  Earth, [Ref,  1] 

The  local  phenomenon  is  comprised  of  those  components  that  are 

not  coherent  over  extended  distances  and  are  assumed  to  be 

products  of  the  local  environn'ient  of  the  detector.  A  set  of 

parallel  detectors  separated  by  some  distance  can  be 

coherently  combined  to  remove  the  global  portion  of  the 

background  noise  and  thus  lower  the  noise  floor. 


B.  OBJECTIVE 


The  Applied  Physics  Lab  (AFL;  at  Johns  Hopkins  university 
has  for  years  studied  and  made  measurements  of  the  background 
geoelectric  noise  on  widely  separated  pairs  of  electrodes. 
The  purpose  of  this  thesis  is  to  describe  techniques  of 
adaptive  filtering  that  the  author  applied  to  data  collected 
by  APL,  in  order  to  determine  the  amount  of  noise 
cancellation.  The  techniques  are  implemented  using  software 
developed  during  an  extended  visit  to  the  Applied  Physics  Lab 
by  the  author.  In  a  related  effort,  software  will  also  be 
pre.sented  implementing  appropriate  detection  m.ethodoiow  for 
the  detection  of  narrowband  sign.ais  embedded  in  the  noise. 


C.  ORGANIZATION 

The  thesis  begins  with 
constituents  of  this  ELF  noise 
in  modeling  of  the  spectrum 
background  infonration  . 


3  brief  discussion  of  the 
to  demonstrate  the  difficuity 
and  to  p rev i  d e  s uppo r  1 1 .n g 

tne  adaptive 


The  basic  concents 
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filtering  v/ill  be  presented  as  a  lead  in  to  the  adaptive 
algorithm  to  be  used  to  cancel  the  noise.  These  introductory 
sections  will  be  follov;sd  by  a  discussion  of  stability  and 
convergence  concerns  for  the  adaptive  algorithm  in  this  type 
of  data  environment .  The  detection  discussion  is  then 
presented  as  a  continuation  of  the  processing  required  for  the 
detection  of  narrowband  signals.  The  adaptive  algorithm  is 
then  used  in  the  detector  scheme  to  demonstrate  performance  of 
the  developed  softv/are.  A  final  section  v;ill  suironarize  the 
work  performed,  draw  relevant  conclusrons  abcuc  the  work  and 
propose  possible  future  studies  that  could  be  derived  from 


this  work. 


II.  THE  ELF  BACKGROUND 


Many  studies  have  been  conducted  in  an  attempt  to  develop 
an  understanding  of  what  actually  constitutes  the  natural  and 
man-made  background  noise  spectrum  in  the  one  to  sixty  hertz 
range.  The  sources  of  the  energy  in  this  frequency  range  can 
be  primarily  attributed  to  atmospheric  disturbances  and  man- 
made  power  transfer  systems.  A  major  problem  with  these 
sources  of.  noise  is  chat  they  are  highly  irregular  in  both 
magnitude  and  duration  of  the  noise  signal.  Attempts  to 
classify  and  assign  specific  characteristics  to  these  sources 
have  proven  difficult.  A  brief  history  of  the  research 
conducted  in  these  areas  is  provided  in  the  following  two 
subsections . 


A.  ATMOSPHERIC  NOISE 

The  source  of  the  majority  of  noise  in  the  one  to  sixty 
hertz  frequency  range  is  lightning  discharges.  The  lower  ELF 
range  discussed  here  also  contains  what  is  called  the  Schumann 
Range.  The  Sch’umann  range  takes  its  name  fromi  the  early 
theoretical  work  of  the  German  scientist  W.  0.  Sch'um.ann . 
Schumann  theorized  that  the  cavity  created  by  the  surface  of 
the  Earth  and  the  Ionosphere  ha.t  naturally  occurring  resonant 
frequencies.  Ke  took  this  resec-rch  further  to  di.ri-.e  a  set  of 
harmonic  resonant  equations  for  cai-ulating  these  frequencies. 


The  v;ork  of  Schumann  has  been  carried  further  with  the 


actual  measurement  of  the  Schumann  resonance  frequencies  in 
field  experiments.  In  1959  and  19G0,  K.  L.  von  Konig 
substantiated  the  presence  of  Schumann  resonances  by  observing 
the  waveforms  in  the  output  of  a  narrov/  band  amplifier.  The 
experimental  presence  of  Schumann  resonances  created  a  large 
quantity  of  literature  on  the  subject  as  scientists  attempted 
to  understand  the  phenomenon.  Once  a  foundation  of  the 
propagation  properties  of  the  Ear th- Ionosphere  cavity  was 
estab’jished,  tlie  search  for  the  source  of  the  excicatior. 


began .  [Kef .  2 J 

The  most  prominent  excitation  source  of  the  Earth- 
Ionosphere  cavity  is  the  cloud  to  ground  lightning  stroke. 
The  number  of  similar  strokes  present  in  any  one  flash 
observed  during  a  thunderstorm  is  a  random  variable  that 
normally  ranges  from  twc  to  twelve.  Since  the  number  of 
strokes  is  generally  too  high  to  n 
a  more  applicable  measure 
density [Hef.  3].  The  frequency 
resolved  in  sufficient  detail  to 
0,  of  the  eurth-ionospherc  ca\  j 
Wagner  in  i960.  The  0  of  a  cavi 
the  bandwidth  of  tt  resonator  [: 
nOiSe,  as  theorized  by  K.  Raemer  in 
the  earth-ionosphere  cavity  to  elec 
in  thunderstorms  all  over  the  wcr 
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model  the  ionosphere  in  an  attempt  to  reproduce  mathematically 
the  measured  response  of  Balser  and  wagner .  Although  Raemer's 
attempt  fell  short  of  its  goal,  it  began  the  process  of 
refinement  of  a  working  model  of  the  ionosphere  that  continues 
today . [Ref .  2 ] 

The  early  work  performed  by  these  scientists  brought  out 
the  difficulty  in  attempting  to  model  this  naturally  occurring 
phenomenon  due  to  its  randomness.  Parameters  must  be 
estimated  by  some  method  in  order  to  model  this  random 
process.  The  nichly  variable  nature  cl  the -e  paran;eters  makes 
the  task  of  modeling  the  noise  nearly  impossible.  Examples  of 
these  varying  parameters  are  diurnal  variations  in  the 
ionosphere,  twenty-four  hour  variations  in  the  ionosphere, 
seasonal  variations  in  the  locations  of  thunderstorm  regions, 
and  the  randomness  of  individual  lightning  strikes.  These  are 
only  a  fev/  examples  of  the  parameters  that  must  be  considered 
in  order  to  model  the  atmospheric  noise. 

3.  MAN-MADE  NOISE 

The  presence  of  atmospheric  ncise  is  not  the  only  concern 
that  iTiUSt  be  addressed  when  discussing  background  noise  in  the 
one  to  sixty  hertz  range.  Many  components  in  today's  world 
emit  unshielded  interference  that  falls  in  this  range.  These 
components  are  the  result  cf  eitiier  faulty  equipment  or  poor 


design.  of  this  nOiS#  SOUiCC©  &£■$  pOWfit  lxn6s  thdt 

ere  used  to  distribute  alternating  current  through  populated 
areas . 

Large  ferric  objects  can  create  fields  around  them  that 
can  be  detected.  The  strength  of  these  fields  can  be  v;ell 
below  the  previously  mentioned  sources  of  noise  but  can  be 
discernable  in  certain  situations.  Motion  of  these  large 
magnetic  objects  in  the  Earth's  natural  magnetic  field  is  the 
source  of  a  small  field. 

Mot  i or,  of  t  lie  detector  can  create  an  increase  in  the 
c' .i  j k  ; I  'i  r'.cise.  Any  siicht  mot  ion  of  the  detector  can 
create  misleading  information  and  in  most  cases,  an  increase 
in  the  noise  signal .  The  background  ELF  spectrum  can  be 
thought  of  as  the  summation  of  three  orthoganai  vectors. 
Relative  motion  between  these  three  vectors  and  the  detector 
creates  a  false  signal  that  can  raise  the  noise  floor.  This 
type  of  noise  can  be  compensated  for  by  the  use  of  motion 
sensor . f  Kef .  5 ] 

wt;',ei.  sources  of  man-rnade  noise  can  be  found  in  the 
detector  equipment  itself.  Improper  shielding  of  cables  and 
electronic  devices  used  in  the  design  of  the  detector  and  its 
processing  componer.tr.  can  lead  to  elevated  noise  floors. 
Unlike  naturally  occurrir.c  atmospheric  noise,  these  sources 
can  be  traced  and  eliminated.  The  removal  of  this  noise  can 
come  in  two  form's,  elimination  of  the  source  or  isolation  cf 


/ 


llx.  ADAPTIVE  FILTERS 


A.  GENERAL  TERMINOLOGY 

The  ELF  background  noise  has  been  described  as  a  highly 
impulsive,  non-stationary  random  process.  Some  of  the  causes 
c£  these  variations  have  been  briefly  described  above.  The 
use  of  traditional  filter  designs  may  not  be  the  optimum 
method  £or  reduction  of  the  noise  floor  for  reception  of  low 
strength  signals  in  the  one  tc  fifty  hertz  range.  Since  their 
introduction,  the  use  ci  adaptive  filters  has  shown  promise  in 
comhating  exactly  this  type  of  scenario.  To  understand  the 
application  of  adaptive  filters  in  this  case,  a  brief 
discussion  of  relevant  terms  must  be  conducted. 

An  adaptiv  '  filter  is  by  definition  a  filter  that  has  the 
capability  to  adjust,  by  itself,  a  set  of  design  parameters 
that  are  based  on  estimated  statistical  characteristics  of 
the  signal  to  be  filtered.  This  process  of  adjusting  the 
filter  to  the  present  situation  alleviates  the  need  tc  have  a 
prion  knowledge  of  the  relevant  signal  characteristics.  The 
previous  chapter  described  how  difficult  it  is  to  derive  a  set 
of  statistical  characteristics  for  the  ELF  noise  environment . 
If  a  priori  knowlec“e  of  the  relevant  signal  characteristics 
were  known,  an  cptimum  filter  could  be  designed.  The  goal  of 
adaptive  filtering  is  to  estimate  these  statistical  parameters 
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and  to  refine  the  estimation  through  the  use  of  an  algorithm. 
If  an  appropriate  algorithm  can  be  found,  then  after  a 
sufficient  numijer  of  iterations,  the  adaptive  filter  should 
converge  to  the  optimum  filter.  The  use  of  adaptive  filters 
requires  the  designer  to  look  for  the  optimum,  algorithm  for 
estimating  the  relevant  signal  characteristics  instead  of 
attempting  to  discover  the  actual  signal 
characteristics . [Ref ,  6] 

Adaptive  filters  are  divided  into  tv.'o  basic  categories, 
open-loop  and  closed-loop.  An  open-loop  configuration  is 
usually  a  tv,-o  stage  process.  The  first  stage  is  used  to 
"learn"  the  statistics  of  the  relevant  signal.  The  results  of 
the  first  stage  are  then  used  in  the  second  stage  to  compute 
the  filter  parameters  required  using  a  nonrecursive  algorithm,. 
In  contrast,  the  closed-loop  configuration  uses  only  one  stage 
to  develop  the  filter  parameters;  the  relevant  signal 
characteristics  are  net  explicitly  estimated  but  are  inferred 
from  a  recursive  algorithm  that  updates  the  filter  parameters 
directly  as  each  new  data  point  is  obtained.  The  adaptive 


filter  gains  additional  knowledge  of 


signal 


characteristics  during  each  iteration.  The  gain  in  knowledge 
results  in  an  improvement  of  the  filter  parameters  v.’hic:;  are 
adjusted  and  their  performance  is  used  a.s  an  input  to  the  next 


iteration  of  the  alaorithm..  The 


•ive  filter  discussed  jn 


this  thesis  fails  in  the  closed  loop  category.  The  advantage 


of  closed-loop  over  open-loop  is  r.hat  the  requirement  for  oniy 


one 


steg©  generally  corresponds  to  a 


i.ess  expensive 


configuration  based  on  its  simpler  implementation . [Ref  6] 

The  description  of  the  filters  that  will  be  analyzed  in 
this  study  carries  teoiinology  that  is  pertinent  to  all 
filters,  not  just  adaptive  filters.  Filters  are  implemented 
in  either  continuous-time  or  discrete-time  which  defines  the 
form  of  the  relationship  between  the  input  and  output  of  the 
filter.  The  filters  that  will  be  studied  are  of  the  discrete¬ 
time  version  v;hich  means  that  the  filter  may  be  described  by 
a  difference  ecuaticn.  T!ie  structure  cf  the  filters  tc  be 
studied  are  referred  to  as  f  in  i  t  e- irr.pih  se  response  {FI?.;, 
tapped-deiay-line,  or  transversal  filters.  This  description 
means  t.hat  the  filter’s  output  relies  only  on  the  past  and 
present  values  of  the  input.  The  advantage  of  a  FIR  filter  is 
its  inherent  stability. [Ref  6’ 

The  algorithm  that  is  used  to  update  the  filter  parameters 


generally  na.med  after  some 


its  operating 


characteristics.  The  algorithm  used  here  is  called 


:ne 


Sequential  Regression  Aigorithiri  or  SER[Ref.  7j.  This 
algorithmi  uses  an  iterative  approach  to  approximating  the 
weighting  factors  required  for  a  linear  combiner.  Figure  1 
shows  a  basic  diagram  of  a  linear  combiner. 

Th-is  system,  uses  simultaneous  sam.ples  of  the  input  signal 
at  time  index  k,  to  produce  an  output  signal,  y.. ,  by  a  linear 
operation  with  weighting  factors,  w.  A  second  way  to 
physically  interpret  the  linear  combiner  is  tc  look  at  it  in 


Figure  1  Basic  Linear  Combiner 


the  form  of  a  transversal  filter.  Figure  2  shows  how  this 
system  is  designed:  corresponds  to  the  last  L  samples  of 
X  at  time  index  k.  The  weighting  factors  are  dual  indexed  on 


the  time  delay,  1,  followed  by  the  time  index,  k.  The  weight 
factor,  W2y,  corresponds  to  the  weight  factor  for  the  sample 
taken  two  samples  ago  from  the  present  sample  time,  k. 


B.  THEORY  AND  DESIGN 


1.  Configuration 

The  filter  to  be  developed  here  is  a  coirbination  of 
the  previous  two  figures.  The  transversal  filter  is  applied 
to  a  set  of  multiple  sensor  outputs  to  form  an  input  vector 
that  is  N  X  L  data  points  long,  v/nere  L  is  the  numbser  of 
separate  inputs  or  sources,  and  N  is  the  length  of  the 
transversal  filter.  There  are  two  reasons  for  using  a  filter 
design  of  this  nature.  The  first  is  the  ability  to  use  more 
than  one  source  of  input,  aliov.-ing  for  a  better  spectral 
sampling  of  ir'.c  background  EbF  encigi'.  The  second  is  the 
ability  to  use  the  previous  N  samples  thereby  smoothing  the 
output  since  it  is  not  dependent  on  instantaneous  values  v;hich 
may  varv'  Greatly  from  sample  to  sample.  Beth  of  these  ideas 
are  beneficial  to  the  overall  performance  of  the  filtering 
process . 

The  output  of  the  filter,  y„,  can  be  expressed  in 
terma  cf  the  input  mar.rix  and  the  weighting  factors  m.atrix. 
To  do  this,  some  defining  equation  m.ust  be  expressed  for  the 
input  matrix ,  X,  , 


^ 

and  v;eight  factors  matrix,  W, , 


••• 


,vhere  each  sensor's  inrut  vector  is 


^  [x,  (if)  x^{k"l)  Xi(if-2)  ...  x^{k~  (N-l))]  (3) 

and  weighting  factor  vector  is 

^ik  ~  ^i{k-Z)  •  •  •  (iif- (W-i) )  5 

The  first  subscript  of  Equations  3  and  4,  labeled  1,  2,  3, 
...,  L,  indicates  the  source  of  the  sample.  Tlie  second 
subscript,  k,  indicates  the  time  of  the  sample.  A  total  of  N 
sample  values  from  each  sensor  are  used  for  each  calculation 
of  y-, .  A  physical  model  for  the  filter  to  be  studied  is  shown 
in  Figure  3.  The  scalar  output  value,  y,  in  Figure  3,  is 
defined  as: 

(5) 

These  basic  equations  will  become  the  starting  point 
for  the  development  of  an  operating  filter  system  for  the 
removal  of  the  ELF  background  noise.  It  is  preferable  that 
the  weighting  factors,  W,.,  cause  the  output,  y>,  to  be  an 


estimiate  of 

some  desired  signal,  d. 

The 

concept  of  the 

filter  is 

to  use  the  weighting 

factors 

to  model  the 

di f  f erences 

between  signals  originati 

ng  from 

the  same  source 

but  sensed  at  various  locations.  The 

inputs , 

X,  in  Figure  3 , 

represent  samples  of  a  signal  taken  at  different  locations. 
The  desired  signal,  d. ,  is  also  a  sample  of  the  same 
background  source  signal.  The  difference  between  the  desired 
signal  and  the  output  of  the  filter  is  expressed  as  the  error 
in  the  estimate,  €,. 
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(6) 


In  the  system  designed  here,  the  estimate  error  is  the 
output  sought  for  analysis.  If  the  weight  factors  are  set 
properly,  they  will  theoretically  remove  all  signals  that  are 
coherent  betv;een  the  different  sample  locations.  The  goal  of 


such  a  system  is  to  remove  the  correlated  background  noise 
from  a  single  sensor  thus  enhancing  signals  that  are 
particular  to  that  sensor  alone.  If  the  estimated  value  of  y^, 
can  be  made  to  represent  the  ELF  background  noise,  then  the 
estimate  error,  can  be  assumed  to  be  due  only  to  the 

signal  present  in  the  desired  signal.  The  resulting  output 
from  the  noise  canceler  has  a  much  lower  noise  floor  due  to 
the  removal  of  the  correlated  noise. 

2 .  The  Mean  Square  Error 

The  filtering  process  must-  have  a  means  by  which  an 
optimal  solution  can  be  defined.  Trie  m.ean-square  error  will 
be  used  to  develop  a  defining  equation,  called  the  performance 
surface,  for  the  optimal  solution.  Combining  Equations  3  and 
4,  a  definition  of  estimated  error  in  tenr;S  of  the  input 
vector,  the  desired  signal,  and  the  -weighting  factors  is 
obtained.  The  expected  value  of  the  squared  estimated  error, 
£(€•,)",  assuming  that  €,,  d.,  and  X,  are  statistically 

stationary  is  given  as: 

E[ei]  =  W  (7) 

The  variables  x,.  and  d,  are  a5sum;ed  to  be  dependent,  and  in 
this  case  both  contain  the  same  ELF  background.  The  eq-aation 
for  the  mean  square  error,  MSE,  associated  v.’ith  the  estimate 
is  thus  defined  in  terms  of  the  auto-corr elat icn  function  of 
the  input  vector  and  the  cross-correlation  fun  c  t  i  o  ri  t)  fc*  V.'  0  Q  n 
the  input  vector  and  the  desired  signal. 


The  autO“Correlation  and 


the  crosS“Corrtlation 


functions  v^ill  be  represented  by  R  and  P  respectively,  The 
elements  of  these  tv;o  matrices  are  constant  values  if  the 
input  vector  x,  and  d,  are  stationary  random  variables.  A 
mathematical  expression  for  the  auto-correlation  matrix,  R, 
is : 

’ 

... 

R  =  E  [X.xD  =E  ■  ■  ‘ 

huXn  ■■■  Xit 

and  the  cross -correlation  n-iatrix,  P,  is; 

P  =  E[d,X,]  =E[dp^,,  ...  d^^]^  (9) 

Using  these  two  definitions  in  Equation  7  results  in  the 
follov.’ing  definition  of  the  mean  square  error  in  the  k-" 
estimate  of  y: 


MSE-J^^  =  E[rf]  =E[. 
The  optimum  solution  for  the 
factors  IS  the  condition  where 
minimum . 

3 .  Minimizing  The  Error 

The  minimum  mean  square 
for  by  taking  the  gradient  of 
vjeighting  factors  and  setting  : 


f|]  +  W^RM  -  2P^W  (10) 

selection  of  the  weighting 
tr.e  mean  square  error  is  at  a 

error  condition  can  be  solved 
the  M.SE  with  respect  to  the 
.t  equal  to  zero.  The  use  of 
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i*' 


the  minirnuni  value  of  the  MSE  as  the  optimal  value  is  based  on 
the  assumption  that  the  weighting  factors  are  only  designed  to 
model  that  portion  of  the  input  vector  that  is  correlated  v/ith 
the  desired  signal.  If  the  error  is  at  the  minimumi  value, 
then  the  filter  is  removing  the  maximum  amount  of  correlated 
data  between  the  input  and  the  desired  .signals.  The  system 
being  evaluated  is  based  on  the  ass'omption  that  only  the 
desired  signal  contains  both  noise  and  a  sign.al  of  interest 
while  the  input  vector  is  composed  of  a  pure  noise  signal 
alone.  Under  this  acsurr.pt  ion ,  the  only  correlated  data 
between  the  two  input  data  sanples  is  the  noise.  The  desire 
is  to  remove  as  m.uch  noise  as  possible  and  thus  to  remove  the 
maximum  amount  of  correlated  data  between  the  input  vector  and 
the  de.sired  signal.  It  is  therefore  assumed  that  the  filter 
is  v/orking  at  its  optimum  value  when  the  weighting  factors 
meet  the  criteria  of  minimizing  the  mean  sguare  error  of  the 
output . 

4.  The  Optirmom  Weighting  Vector 

The  optim.'UiT'  weight  factor  will  be  represented  by  the 
syniocl  W  .  The  two  conaiticns  that  m.ust  be  met  in  r^rder  for 
the  value  of  J  to  be  minimized  are: 

^65  ® 

and 

.  e 


dwj^dwj 


These  equations  state  that  vjhen  the  v/eight  vector  is  at  an 
optimum  the  perfonuance  function  is  at  a  local  minimum  since 
the  slope  is  nero  and  the  curvature  is  upv/ard.  The  gradient 
of  the  performance  function  is  evaluated  to  be: 

V  =  %  J=V_£[d;f]  -2%(W^P)  +y„(W^RW) 

(13) 

=  -2P  +  2RW 

Setrina  the  gradient  to  cero  and  solving  for  the  optimum 
weighting  factors  results  in: 


to  -  B-l 


R-*P 


This  equation  is  sometimes  referred  to  as  the  Weiner-Hopf 
equation  where  the  weight  factor  is  referred  to  as  the  Weiner 
weight  vector.  The  input  autocorrelation  matrix,  R,  can  be 
inverted  if  it  is  "positive  definite,"  i.e.,  if; 


V^RV  >  0 


where  V  is  any  non-zero  vector.  When  the  matrix  R  is 


)sit  ive 


definite"  then  the  optimum  weight  factors  can  be  solved  for 
directly  if  the  complete  statistical  properties  of  the  auto¬ 
correlation  and  cross-correlation  functions  are  known.  The 
statistical  properties  are  not  known  in  this  case  so  further 
development  is  required  to  obtain  an  estimate  of  the  ootimum 


wGignt  factors. 
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Tha  second  condition  that  ensures  that  the  error 
estimate  is  at  a  minimuiTi  is  that  the  partial  derivative  of  the 
gradient  v;ith  respect  to  the  v;eighting  factors  is  greater  than 
zero.  To  prove  this  the  we  take  the  partial  derivative  of 
Equation  3  with  respect  to  W  ; 

=  2R  (16) 

dWjdWj- 

This  result  again  shov/s  that  the  input  autocorrelation  matrix 
nusc  be  "positive  definite"  in  order  to  solve  for  the  optimum 
vveichtinq  factors. 

The  perform.ance  surface  can  now  be  expressed  in  terms 
of  the  minimum  error  possible,  the  weighting  functions , and  the 
autocorrelation  of  the  input  vector.  The  minimum  mean  square 
error  possible  can  be  found  for  the  case  when  the  weighting 
vector  is  at  its  optimum  value.  The  optimum,  v/o.; -iliting  factor 
is  used  to  solve  for  the  minimum  error  by  substituting  W  into 
the  m.inim.um  error  equation. 

=  E[d^]  + 

=  r  [R-^P]  “RR-P-2P"R-1P 

=  E[d^]  -  P^R-'P 
=  E\dl]  -  P^W^ 


The 

mdnimum  error 

can 

now  be  u 

sed  as  the  base  line  by  w 

hich 

al  1 

other  errors 

are 

found . 

The  error  at  any  time  is 

the 

min; 

Lmum  error  wi 

r.h 

the  addi 

ticn  of  uncertainty  in 

the 

components  used  to  derive  the  minirr.'um  error.  Thus: 
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(18) 


The  expression  can  be  sho\-ni  to  equal  that  of  Equation  10  by 
the  following  proof.  Using  the  fact  that  (A-B)"  in  general 
can  oe  expressed  as 

J  -  [WO]  [W°]^’'RW  (19) 

since  all  terms  are  scalars,  the  transpose  is  equal  to  the 
scalar  val'.ie,  thus  the  last  tv;o  terms  are  equal.  Substituting 
for  the  ininimurr  error  and  the  optimuiri  weighting  vector: 

J  -  E[aJ]  -P^fe‘‘P+  [H"]  "RW"'+W^RW-2W^RW" 

=  E[d-^]  -P^-^R  +  P^-*RR"^P+1f^RW-2W^RR-'P 
=  E{d'h  +W^RW-2W^P 
=  i,’[d;?]  +W^RK-2P^ 

This  result  validates  Equation  18  by  showing  that  it 
equivalent  to  Equation  10.  An  important  note  about  this 
expressioii  of  the  perfcrm.ance  surface  is  that  it  is  quadratic 
with  respect  to  the  v.^eighting  vector. 

5.  Iterative  Calculation  of  V/, 

The  optimum  weight  factors  can  be  arrived  at 
iteratively  by  use  of  a  gradient  search  of  the  per f 'urmance 
function.  ’What  this  mo-ana  is  that  after  each  estimate  of  the 
optimum  v.-eight  vector,  W. ,  the  slope  or  gradient  of  the 
pe; formance  function,  in  the  direction  of  the  optimum 
weigiiting  vector,  is  used  to  derive  the  next  Vv'eighting  vector, 
W. , .  .  This  method  ensures  that  each  estimation  of  the 


ootimurr. 


weight  factors  is  at  least  as  close  or  closer  to  the  optirnuiti 
value  than  the  previous  estimate.  Multiplying  Equation  13  on 
the  left  by  H  'R''  an-5  combining  the  result  with  Equation  14, 
we  obtain! 

\)o  ^  (21) 

2 

To  convert  this  into  an  iterative  process,  the  optimum, 
solution  is  considered  to  be  the  value  of  the  v;eighting  vector 
for  the  next  iteration. 

Wa-i  <22) 


This  equation  is  referred  to  as  Che  Newton  Method  algorithm, 
for  determining  the  weigiiting  vector. 

The  algorithm  takes  only  one  iteration  to  arrive 
theoretically  at  the  optimum  solution  for  the  values  of  the 
weighting  vector.  This  algorithm,  would  be  the  ideal  method  to 
u.se  if  exact  solutions  to  the  gradient  term,  and  the  inverse 


correlation  of  the  input  vector  were  kr.O'v\m  exactly.  This  is 
not  generally  the  case  and  an  estim.ate  for  the  gradient  must 
be  used.  To  compensate  for  the  lack  of  knowledge  of  the 


gradient  at  iteration  k,  a  constant  u,  will  be  used  instead  of 


the  factor  hi  in  Equation  22.  The  requirements  on  the  range  of 
values  shat  u  can  take  on  are: 


0  <  <  1 


(23) 


The  effect  of  u  on  the  solution  at  each  iteration  is  that  it 
affects  the  rate  cf  convergence  of  the  algorithm.  A 
discussion  of  the  attributes  and  further  limitations  on  the 
value  of  u  will  be  covered  later.  A  general  form  of  the 
Newton's  method  of  gradient  search  is  therefore: 


(24) 


This  expression  is  ideal  under  the  following  three  conditions: 


i  .  p  =  id 

2.  Exact  knov.'ledge  of  the  gradient  vector, 

3.  Exact  knowledge  of  the  (unchanging)  matrix  R‘' . 

These  conditions,  unfortunately,  are  not  normally 
attainable.  The  value  c£  u  is  selected  between  C  and  H  so 
that  the  algorit'nm,  is  overdamped  and  can  accommodate 
fluctuations  in  the  matrix  R.  The  compensation  for  this 
.selection  of  u  is  that  the  alyoritham  takes  longer  to  obtain 
its  optimu.m  solution.  A  second  modification  to  the  value  of 
u  is  to  base  its  value  or.  the  average  eigenvalue,  A,,,,,  so  that 
u  is  replaced  by  uA_,.. . 


to  the 
at  taina 


manner . 


The  second  condition  pertaining  to  the  exact  solution 
gradient  of  the  error  surface  is,  in  general,  never 
ble.  The  gradient  must  therefore  be  estimated  in  som.e 
The  estimate  of  the  gradient  that  will  be  used  is 
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based  on  the  square  of  the  present  value  of  the  error  signal 
and  not  the  true  ensetrhle  average  as  it  has  been  up  to  now. 


From  the  definition  of  the  error  signal: 


(25) 


(26) 


which  leads  to  the  following  estimate  of  the  gradient: 

-  -2e»I,  (27) 

With  these  two  relaxations  from  the  three  ideal  conditions, 
the  algorithm  is  now  in  a  form.  Icnown  as  the  Kewton/LMS 
algorithm.fRef  7]  : 

"k-l  (28) 

This  is  an  optimum,  iterative  algorithm  that  is  only  com.plete 
if  full  knowledge  of  the  input  correlation  function  is  known. 
This  is  generally  not  the  case  and  thus  a  final  modification 
m.ust  be  done  to  produce  a  working  estimate  to  this  algorithmi. 
The  final  m.odif ication  is  an  iterative  estimate  of  the  inverse 
of  the  input  correlation  matrix,  R'^,  and  is  called  the 
Sequential  Regression  (SER)  algorithm. 


IV.  THE  SEQUENTIAL  REGRESSION  (SER)  ALGORITHM 


A.  DEVELOPMENT 

The  Sequential  Regression  (SER)  algorithm  uses  the 
LMS/Newton  algorithm,  as  its  basis  in  conjunction  with  a  method 
for  iteratively  estimating  the  value  of  the  inverse  input 
auto-correlation  matrix  to  produce  the  weighting  vector  for 
the  filter  coefficients.  The  development  of  the  algorit'nm. 
comes  from.  Adaptive  Sigr.al  Prccessir.c  by  Widrow  and 
Stearns  [Ref  7]  .  The  core  of  the  algorithm,  is  the  use  of  an 
estim,ate  for  the  inverse  of  the  input  auto-correlation  matrix 
which  reduces  com,putational  load  with  ir.ini.mal  estimation 
error . 

The  input  auto-correlation  matrix  is 


R=E[X;.Sj]  (29) 

v;here  the  subscript  k  goes  over  the  entire  lencth  of  the 
random  process  X  .  R  will  be  estimiated  iteratively  using  on'y 
a  finite  or  truncated  sample  of  the  random  process  X..  The 
estimate  of  the  input  auto-correlation  matrix  is  given  by: 

a.-  =  I3U 

A  +  *  2»C 

This  estimate  is  unbiased  when  the  input  variable,  X; ,  is  a 
stationary'  random  process.  When  X.  is  not  stationary,  it  can 
be  shO'//n  that  the  estimate  is  not  very  accurate.  The  estimate 
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of  R,,  at  each  value  of  k  is  equally  influenced  by  the 
previous/  k-l,  estimates  of  the  input  auto-correlation  matrix. 


In  order  that  the  most  recent  estimates  of  Rv,  have 
influence  over  the  new  estimate,  a  "fading"  memory  term,, 
introduced  which  reduces  the  significance  of  the 
estimates.  The  fading  memory  expression  of  the  input 
correlation  matrix  is  denoted  by  where: 

k 


Qk  = 


a 


k-l 


more 
«/  is 
past 
auto- 


(31) 


The  value  of  a  is  chosen,  as  a  rule  of  thumb,  such  that  the 
half  life  of  the  exponential  function  is  equal  to  the  length 
of  stationarity  of  X;[Ref  7]  .  This  rule  of  th’orh  leads  to  the 
following  statements  about  the  value  of  a. 


(32) 


0  <  a  <  1 


(33) 


"he  value  of  the  fading  memory'  factor  ever  k  iterations  is: 


-a 


k*l 


1-a 


(34) 


and  thus  the  estim.ate  of  the  input  auto-correlation  m.atrix 
becomes : 


*  1-a''*""'' 


T  a*-^ 


z,x 


1-^1 


(35! 


2  b 


This  estimate  is  exact  for  the  condition  where  is  constant. 
If  the  inp'jt  variable,  X^,  is  stationary,  a  is  one,  and  if  the 
limit  of  Equation  35  is  taken  as  a  approaches  one,  the  result 
is  in  agreement  with  Equation  30. 

The  derivation  of  the  SER  algorithur,  begins  with  the 
estimate  of  the  input  autC“Ccrrelation  matrix.  Using  this 
estimate  in  the  equation  for  the  optimum,  weighting  factors 
results  in 


(36) 

The  definition  for  tne  estim.ate  of  the  cross-correlation 
vector,  Py, ,  is,  similarly  to  Equation  35, 


r  a 


k~ 


(37) 


Using  Equations  35  and  37  in  36  gives: 


(38) 


The  SER  aigorithmi  atte:rjt.s  to  estimate  the  next  value  of 
the  v.’eighuing  vector,  W,.,  ,  based  or.  the  present  values  of  R, 
and  P.. .  To  do  this.  Equation  38  becomes; 
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(39) 


i“0 

-  aQ^f.j^Wjj  +  djjXjj 

The  last  line  c£  Eq^aation  3  9  is  the  result  of  the  following 
relationship : 

Q,  =  t,xl.  (40) 

The  value  of  d..  can  be  replaced  by  recalling  Equation  5  and  6 
which  state  that 


=  H;,  +  (43 ) 

The  e>:pression  for  W^.-  in  Equation  43  is  similar  to  the 
LMS/Newton  method  derived  earlier.  With  this  in  m.ind,  an 
approximation  for  the  LKS/Newton  algorithm,  can  be  made  using 
the  above  derivation.  The  SER  algorithm,  used  as  an 
approximation  of  the  LMS/Newton  m;ethod  is  given  as: 


Ar.  iterative  method  of  calculating  Q,  "  is  nov,-  available  and  is 
only  based  on  its  previous  value,  Q..  and  the  present  value 
of  the  input,  X...  The  value  in  parentheses  is  used  three 
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it  will  be  calculated  separately  and  represented  by 


times , 

s,. 

The  SER  algorithm  is  now  completely  derived  as  an 
iterative  method  of  approximating  the  LNS/Newton  method  of 
adaptive  filtering.  The  following  summary  show.s  the  process 
by  which  each  nev;  value  of  the  weighting  factor  vector  is 
calculat  ed . 


^  jw  2 stationarity  of  X) 

Y  *  a+xjs 


0  <  ^  <  ^ 
^max 


(49) 


The  values  that  must  be  carried  forward  from,  one  iteration  to 
the  next  are  Q.  -'h  and  W, .  These  represent  the  previous  value 


for  the 

e  E  t i ma  r 

e  o 

f  the  inverse 

input  a':^ 

to-ccrrel 

at  ion  m.atrix 

and  the 

present  es 

timate  for  th 

G  v;eight 

ing  facto 

r  vector 

The 

output 

of 

the  algorithm': 

is  t he 

srror,  c,, 

between 

'.:hG 

present 

value 

of 

the  desired 

signal 

and  esti 

ma  t  e  of 

desired 

signal 

(50) 


OUTPUT  ^  ^ 

This  value  is  saved  after  each  iteration  in  a  vector 
representing  the  uncorrelated  portion  of  the  desired  signal 
and  noise  vector  of  the  primary  channel  with  the  noise  only 
vectors  of  the  reference  channels  at  the  time  index  k, 

Xv.  [Ref  7  j 


B.  INITIALIZATION 

The  SER  alcorithur.  assumes  that  the  value  of  the  inverse 
input  auto-correiaticn  matrix  and  the  -weight  factors  vector  is 
known  frcr;.  the  nrevicus  step.  At  some  point  in  the  past,  an 
estimate  for  these  term.s  must  have  been  made  so  that  a 
starting  point  can  be  established.  The  required  accuracy  of 
the  estim.ates  depends  on  the  length  of  time  that  the  algorithm 
will  be  operating  and  the  earliest  time  that  accurate  data  is 
required . 

The  driving  force  behind  the  accuracy  of  the  initial 
estimate  is  the  speed  with  which  the  algorithm,  can  converge  to 
an  optimal  solution.  The  speed  of  convergence  is  dependent  on 
two  eiem.ents  of  the  algorith.rr,.  The  first  element  is  the 
initial  value  of  the  inverse  input  auto-correlation  matrix, 
Q,  The  significance  of  this  initial  element  is  that  the 
further  the  initial  estin;ate  is  from  the  true  value,  the 


rurtner  tne  a^g 


;  c*  ^ 


adac-t  in  order  to  aooroximate  the 


optimum  solution,  The  second  elemient  is  the  step  size, 


which  specifies  how  far  the  aigorithjn  can  progress  towards  the 
optimum  value  in  one  iteration.  The  step  size  will  be 
discussed  in  the  next  section  but  for  now  it  shall  be  assumed 
to  be  a  constant  in  this  discussion  of  the  initial  conditions. 

The  use  of  a  large  constant  multiplied  by  the  identity 
matrix  as  the  initial  value  of  the  has  been  sufficiently 
argued  by  Lee  [Ref.  8]  . 

Qo*  -  q^l  (51) 

As  discussed  by  Lee,  the  use  of  a  large  constant  for  q  is  the 
proper  selection  for  the  case  vinere  the  random  process  is 
Stationary.  Widrov;  and  Stearns  argue  that  this  same  choice 
for  the  initial  value  of  q.-  iR  ^  non-stationary  condition,  is 
also  appropriate [Ref  7]  .  The  point  is  made  that  if  a  priori 
knowledge  of  the  random  process  allows  for  a  better  estimate 
of  the  value  of  Q:'h  then  that  information  should  be  used. 
The  method  of  initialization  for  this  study  is  to  use  Equation 
with  q,  equivalent  to  IOC. 

The  initial  values  used  in  the  weighting  factor  vector  are 
ass'.rntd  to  be  zero.  This  ass'umption  is  made  since  knov;iedce 
of  the  process  before  tim.e  zero  is  unknown.  The  use  of  a  32 
tap  delay  filter  also  allows  for  an  estimate  of  the  initial 
full  weight  factor  vector  but  this  is  done  so  at  the  expense 
of  the  first  32  output  data  points.  This  is  done  by  expanding 
the  numher  of  eiem.ents  in  the  weight  factor  vector  jS  each  of 
the  first  32  points  are  read  into  the  algorithm.. 


aigoi'ithiTt  is  applied  only  to  the  data  chat  is  availabl#; 
therefore (  the  first  time  that  each  element  of  %  is  non-zero 
is  at  time  index  31.  This  corresponds  to  the  time  when  the 
first  32  data  points  are  input  into  the  algorithm.  The  result 
of  this  initialization  is  that  the  first  real  data  point  to  be 
considered  an  output  of  the  SER  method  used  here  is  at  k  s  31. 
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V.  STABILITY  AND  CONVERGENCE 


A.  THE  SINGLE  WEIGHT  CASE' 

The  SER  algorithm  contains  values  that  can  be  set  to 
affect  the  rate  at  which  the  algorithm,  achieves  the  optim.urr, 
solution.  A  simple  case  where  only  one  weight  factor  is 
concerned  will  be  used  to  demonstrate  the  relationships 
involved  in  the  selection  of  tlie  convergence  rate  factors  . 
The  single  weight  iterative  process  using  a  gradient  search 
metiiod  like  the  SER  algorithm  can  be  represented  by: 

w,.,  =  (S2) 

The  expression  for  the  single  weight  gradient  at  tim.e 
index  k  is  given  by: 


dw 


=  2X  {w.^-w°) 


(53) 


The  value  of  the  input  auto-correlat ion  matrix,  r.-  for  the 
single  w'eight  case,  is  replaced  by  its  equivalent  eigenvalue, 
k  =  E[x-]  .  Combining  Equations  52  and  55,  an  iterative  equatior: 
for  the  weight  factors  is  found; 


54 


^k*i  ~  Wjt-2|iX(iirjj.-v°)  (54) 

for  the  single  weight  case.  The  terms  in  Equation  54  can  be 
rearranged  to  combine  like  terms. 

Vj^=w^+{1-2\iX)^{Wq-w^)  (55) 

This  expression  contains  the  geometric  ratio  oi  successive 
terms  which  is  also  used  to  define  the  stability  requirements. 

Stability  is  guaranteed  if  the  absolute  value  of  the 
geometric  ratio  is  less  than  one,  i.e., 

ll-2jiX|<l  (56) 

The  limits  placed  on  the  value  of  the  p  can  be  expressed  in 
terms  of  the  eigenvalue  of  the  input  auto-correlation  matrix 
k  by  rearranging  Equation  56 

-^>|i>0  (57) 

The  optimum  solution  can  be  seen  as  the  solution  of  Equation 
55  if  the  limit  of  as  k  goes  to  infinity.  This  is  because 
the  geometric  ratio  approaches  zero  as  k  approaches  infinity. 
The  rate  at  which  the  algorithm  approaches  the  optimal  value 
is  dependent  on  the  value  of  the  geometric  ratio. 

Figure  4  shows  the  effect  of  varying  values  of  the 
geometric  ratio  on  the  rate  of  convergence.  If  the  value  of 
the  geometric  ratio  can  be  maintained  between  zero  and  one, 
then  the  algorithm  will  always  approach  the  optimal  solution 
from  one  side  and  can  be  considered  as  overdamped.  For  a 
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geometric  ratio  of  zero,  the  algorithm  reaches  the  optimum 
value  in  one  iteration  from  one  side  and  can  be  considered  as 
critically  damped.  If  the  magnitude  of  the  geometric  ratio 
falls  between  zero  and  negative  one  then  the  algorithm  is 


Convergence:  overdamped  0<g<l;  critically  damped  g=0; 
underdamped  -l<g<0 

B.  MULTI -WEIGHT  CONDITION 

The  single-weight  condition  discussed  in  the  previous 
section  can  be  expanded  to  develop  an  understanding  of  the 
multiple  weight  factor  case.  The  same  conditions  that  applied 
to  the  single  weight  case  are  true  in  the  multiple  weight 
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case.  The  difference  is  that  each  weighting  factor  has  an 
associated  eigenvalue  and  thus  the  condition  of  stability  is 
now  dependent  on  a  set  of  simultaneous  equations  that  must 
meet  the  stability  condition. 


(i-2^iA2) 


V 


(58) 


( 


The  subscript  on  X  indicates  the  index  value  corresponding 
a  specific  weight  value.  There  are  a  total  of  L  weighti 
factors  and  thus  there  are  L  individual  eigenvalues,  X. 

The  value  of  u  in  each  of  the  simultaneous  equations 


to 

-g 


is 


held  constant.  Since  u  is  the  same  for  all  equations,  the 
eigenvalue  that  has  the  largest  value  will  iiripose  the  most 
restrictive  conditions  on  the  selection  of  a  value  for  p,  The 
r.ev,»  range  of  values  t.hat  u  can  take  on  is: 


(59) 


The  ccr.dicicr. 
stability  reqv. 
algorithm  will 
The  rate  a 
Figure  4  tc  be 
New  that  the 
eigenvaj. ue  or 


stated  in  Equation  59  gives  the  mult  '-weight 
irements  that  must  be  met  to  ensure  that  the 
converge  o  the  optimal  solution, 
which  the  algorithun  will  converge  was  shov;n  in 
dependent  on  the  value  of  the  geometric  ratio, 
value  of  u  is  a  function  of  the  maximum 
Ine  input  autc-correlat ion  matrix,  R,  it  can  be 


37 


ihown  how  this  affects  the  rate  of  convergence.  By  holding 
the  value  of  u  constant  for  each  weighting  factor,  an 
undesirably  slow  convergence  rate  can  occur  for  widely  varying 
values  of  A,.  As  shown  previously,  the  algorithrr.  converg(!S  in 
one  iteration  when  u  is  egual  to  1/2X  for  the  single  weight 
case.  Using  that  same  criteria,  the  value  of  the  multi-weight 
constant  u  is  l/2Ar.,^.  The  use  of  this  criterion  ensures  that 
the  algorithm  falls  into  the  overdamped  category  since  u  is 
always  less  than  or  equal  to  1/2A,;  for  1=[0,  1,  2,  .  .  .  ,  L;  . 
The  weight  factor,  w  .  with  the  smallest  value  cf  A,  v/ili  have 
the  slowest  convergence  rate  since  its  geometric  ratio  will  be 
farthest  from  zero.  This  is  true  because  the  value  of  2uAi 
will  be  at  its  .minimum,  and  will  cause  the  geometric  ratio  to 
approach  one.  When  this  condition  holds,  the  algorithmi  may 
take  excessively  long  periods  of  tim.e  to  converge  to  an 
optimal  solution. 

C.  NON-STATIONARITY  EFFECTS 

Th.e  discussion  up  to  this  point  in  dealing  with  stability 
and  convergence  ass’om.ed  that  the  process  under  study  was 
stationary.  Under  stationary  conditions,  the  \alu9S  of  the 


eigenvalues 

of 

R  are  assumed  to  be 

constant . 

A  .second  key 

assumption 

to 

the  deveiopmie'it  is 

that  the 

values  of  the 

eigenvalues  are  known.  The  lack  of  knowledge  of  the  input 
auto-correlation  matrix,  R,  requires  that  some  estimation  be 
made  for  the  value  of  y.  The  value  of  u  must  be  adjusted  to 
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ensure  convergence  of  the  algorithm  in  an  environment  where 
the  eigenvalues  of  the  auto-correlation  matrix  are  varying  at 
some  unknown  rat©  from  iteration  to  iteration. 


A  method  for  compensating  for  the  variations  in  the  range 
of  values  that  u  can  assume  is  to  use  an  averaged  eigenvalue, 
for  the  selection  of  the  constant  u  for  each  iteration. 
The  value  of  u  is  now  a  constant  for  a  single  iteration  but  is 
changed  from  iteration  to  iteration  based  on  the  average 
eigenvalue.  Stability  of  the  algorithm  is  guaranteed  if  the 
value  of  the  product,  is  betv.'ee-;  zero  and  one.  It  is 

important  to  note  that  the  expression  uA-av.-  is  now  used  only  to 
replace  the  constant  u  used  previously.  The  idea  is  to  scale 
the  average  eigenvalue  for  that  iteration  by  a  constant,  u, 
and  use  it  as  the  factor  for  controlling  convergence  of  the 
algoritbn. 

The  geo.metric  ratio  which  controls  the  convergence  of  the 
algorithm,  has  now  beer,  altered  to  include  a  term  to  compensate 
for  the  ncr.-stationarity  of  the  input  random  process,  X.  The 
new  expression  for  the  geometric  ratio  is 


|1  -  <1  (60) 

The  inverse  input  auto-ccrrelation  matrix,  R'-,  is  included  in 
the  development  of  the  SER  algorithm.  The  converge.nce 
controlling  term,  is  thus  and  must  be  less  than  the 
maxim.am  eigenvalue  of  R  to  ensure  stability  and  convergence. 
A  second  condition  that  can  be  im.posed  on  the  product  uX.,.,.  is 
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that  its  value  must  always  be  between  zero  and  one.  Using  the 
fact  that  the  spread  between  ma>cimuin  and  minimum  eigenvalues 


can  be  rather  large,  the  product  of  must  be  much  less 

than  one. 

The  ELF  background  noise  processes  have  shown  a  tendency 
to  have  a  wide  range  between  eigenvalues  and  to  slowly  change 
over  time.  Because  of  these  two  conditions,  the  value  of  the 
must  be  adjusted  over  time  to  account  for  these 
variations.  The  chosen  method  is  to  estimate  the  average 
eigenvalue  by  taking  the  average  of  the  square  of  the  input 
vector,  X„.  This  value  is  multiplied  by  the  1/100  of  the 
inverse  of  the  maximum,  of  the  squared  input  vector. 


100  [max  (Xj)  ]  ^ 


(61) 


The  use  of  this  estimate  for  convergence  control  in  the 
algorithuTi  ensures  that  the  requirements  for  stability  are  met 
for  each  iteration.  The  1/100  term  ensures  that  the  maximum, 
value  of  y  is  not  exceeded.  It  was  derived  empirically  with 
the  use  of  an  actual  set  of  EL?  background  measurements  and 
may  require  adjustment  in  a  different  environm.er.t . 

The  choice  of  using  the  present  input  to  the  algorithm,  to 
get  an  estimate  of  the  average  eigenvalue  is  used  for 
computational  ease  with  m.inimal  added  error.  The  use  of  the 
inverse  input  auto-correlation  matrix  should  provide  a  m.cre 
accurate  measure  of  the  eigenvai'ues ;  That,  however,  requires 
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the  inversion  of  the  matrix.  The  prevention  of  matrix 
inversion  is  the  exact  reason  the  SSR  algorithm  is  developed/ 
80  to  obtain  the  eigenvalues  in  this  way  defeats  the  reason 
for  the  algorithm^' s  design,  and  is  therefore  counter¬ 
productive.  The  method  has  proven  adequate  and  reasonably 
accurate  in  empirical  work  up  to  this  point,  but  does  leave 
room  for  improvement  and  further  study. 
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VI .  IMPLEMENTATION 


A.  A  PHYSICAL  DESCRIPTION 

The  Applied  Physics  Lab  has  two  identical  sensor  platforms 
submerged  in  shallow  water.  Each  platform  has  electrode  pairs 
for  detecting  the  geoelectric  ELF  background  noise.  The 
system  has  been  in  operation  for  over  two  years  collecting 
samples  of  the  background  noise  environment. 

Ideally,  the  two  platforms  are  parallel  when  the  piatformiS 
are  placed  on  the  bottom  of  the  bay.  Figure  5  shows  the  ideal 
alignment  of  the  platforms.  The  two  platforms  are  labeled  as 
the  east,  E,  and  west,  W,  platforms.  Each  platform,  has  a  set 
of  orthogonal  electrode  pairs,  labeled  X  and  Y.  These 
electrode  pairs  are  orientated  to  receive  the  tv;o  orthogonal 
source  vectors,  E*  and  E^.,  in  the  horizontal  plane.  The 
strongest  geoelectric  field  strengths  are  in  the  horizontal 
plane  due  to  the  shift  in  polarization  as  a  wave  crosses  tr.e 
water-air  boundary  from  vertical  to  horizontal. 

If  the  two  platforms  are  perfectly  aligned  and  the  noise 
recorded  on  each  are  identical,  the  noise  in  a  specific 
direction,  X  or  Y,  can  be  removed  by  subtraction.  To  cancel 
the  x-direction  noise  simply  subtract  XW  from  XE .  Because  the 
platforms  are  not  perfectly  parallel  and  the  recorded  .noises 


are  not  ident?.cal  though  highly  correlated,  this  method  is  not 
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Figure  5  Illustration  of  Sensor  Set-L’j 


ideal.  The  actual  platforcis  can  be  misaligned  up  to  1 
degrees  from,  perfectly  parallel.  This  is  due  to  physical 
constraints  such  as  bottom  contours  and  mounting  technique. 
The  result  of  such  a  misalignment  is  that  i^oth  the  XW  ar.d  YV; 
sensors  of  the  west  platform  have  components  that  are  coherent 
with  XE  of  the  east  platforrri.  If  the  electrode  pairs  in 
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Figure  are  thought  of  as  vectors,  then  it  is  easily  seen  how 
misalignment  can  cause  coherency  between  XE,  XW,  and  YW. 


A  set  of  real  data  is  traced  through  the  signal  processing 
system  designed  to  accompany  this  sensor  system.  A  calibrated 
source  signal  is  embedded  in  the  signal  XE.  The  path  that  the 
received  signal  follows  can  be  seen  in  Figure  6.  The  output 
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of  aach  sensor  is  sampled  at  128  Hz.  The  data  is  bandpass 
filtered  and  resampled  down  to  32  Hz  for  processing.  The  data 
is  then  passed  through  the  SER  algorithm  where  the  coherent 
portions  of  the  v/est  platform  are  removed  from  XE.  The  output 
of  the  SER  portion  of  the  processor  is  demodulated  and  passed 
through  a  likelihood  ratio  detector  where  a  decision  can  be 
made  as  to  when  the  embedded  signal  is  present. 


B.  INPUTS  TO  THE  SER  ALGORITHM 

The  SER  algorith.m  is  used  to  reduce  tr.e  ELF  background 
noise  detected  in  the  primary  submerged  electric  field  sensor 
by  cancellation  with  the  output  of  the  reference  parallel 
electric  field  sensor.  These  parallel  reference  sensors  are 
also  submerged  and  are,  ideally,  assum.ed  to  be  located  within 
the  same  ELF  background  noise  as  the  primary  channel.  The 
primary'  channel  y  represents  sensor  XE  and  the  reference 
channels,  xl  an.!  x2,  represent  XW  and  Tx',  respectively.  The 
samples  are  the  actual  output  of  the  sensor  system;  used  by 
APL.  The  sam.ples  are  taken  at  the  same  time  without  any  knov.T; 
signals  prese.nt .  The  sampling  rate  of  the  data  input  for  the 
SER  algorithm  is  32  Kz  and  the  total  length  of  data  is  113760 
points.  These  data  samples  correspond  to  almost  an  hour  of 
background  noise  measurements. 

The  filter  is  a  32  point  transversal  filter  that  is 
applied  to  each  channel  simultaneously.  The  length  of  the 
transversal  filter  m.eans  that  the  previous  31  samples  of  the 
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iensor  are  used  with  the  present  sample  value  for  each 
iteration  of  the  algorithm.  The  goal  of  the  algorithm  is  to 
take  the  output  of  the  two  reference  channels,  xl  and  k2,  and 
to  generate  an  estimate  of  the  noise  in  the  primary  channel  y. 
Figure  7  shov/s  what  the  primary  channel's  data  looks  like  in 
real  time.  The  length  of  stationarity ,  used  to  determine  the 
forgetting  factor,  a,  is  assumed  to  be  about  300  data  points 


Figure  7  The  Priniary  Noise  (Y)  and  Calibrated  Source  Signal 

long.  This  value  is  used  based  on  empirical  analysis  of  the 
real  data. 

The  effectiveness  of  the  algorithm  is  analyzed  by  running 
a  sec  of  scenarios  with  a  knov^m:  signal  inserting  in  the  data 
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samples.  The  scenarios  consist  of  manually  inserting  a  32  Hs 
signal  on  a  5  H2  carrier  that  is  8192  data  points  long. 
Figure  7  shov;s  the  components  of  the  primary  channel  signal  in 
a  real  time  representation.  The  primary  channel  is  processed 
by  the  SER  algorithun  while  varying  the  strength  and  location 
of  signals  embedded  in  the  background  noise  y. 


C.  SIGNAL  STRE3JGTH  RESULTS 

The  algorithm,  is  first  run  with  the  location  of  the  signal 
within  the  noise  held  constant.  The  strength  of  the  signal  is 
varied  by  dividing  each  data  point  in  the  signal  vector  by  a 
constant.  A  series  of  four  runs  are  conducted,  each 
succeeding  run  has  half  the  signal  strength  of  the  preceding 
run.  The  initial  run  uses  the  signal  shown  in  Figure  7  with 
the  amplitude  divided  by  500.  The  signal  is  indistinguishable 
prior  to  processing  as  seen  in  Figure  G. 

The  signals  for  each  run  are  centered  at  time  19.2,  32  and 
44.8  itinutes  into  the  noise  sa.mple  y.  There  is  no 
significance  to  the  positions  chosen  except  that  the  signal 
data  streamis  were  not  allowed  to  overlap.  The  objective  of 
this  first  set  of  runs  is  to  demonstrate  how  the  SER  algorithm, 
enhances  the  ab:' iity  to  distinguish  5m.all  signals  buried  in 
the  noise.  The  data  is  first  processed  in  the  SER  algorithm 
to  remicvG  the  coherency  between  the  reference  channels  and  the 
primary  channel.  The  output  is  then  demodulated  leaving  only 
the  in-phase  and  ouadrature  comiDcnents  of  the  non-coherent 
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portion  of  the  primary  channel.  An  identical  reference  set  of 


data  is  demodulated  without  passing  it  through  the  SER 
algorithm  so  that  a  visual  and  quantitative  analysis  of  the 
effects  of  the  SER  algorithm  can  be  made.  The  output  of  the 
demodulator  is  showai  in  Figures  9  and  10.  Note  that  the 
envelope  of  the  signal  is  visible  to  nearly  the  strength  of 
1/2  00  0''',  .  :  the  third  run  made. 

The  output  of  the  demodulator  is  passed  as  input  into  a 
matched  filter  detector  based  on  the  demodulated  signal.  The 
outout  of  this  derector  is  shown  in  Figures  11  and  12  with  the 
SER  and  non-SER  output  for  a  given  signal  strength  shovm 
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together.  Numerically,  the  reak  value  out  of  the  detector  is 
three  to  five  times  greater  when  the  SFR  algorithni  is  used  as 
compared  to  v;hsn  it  is  not  used. 

A  final  numerical  analysis  is  perform, ed  by  using  a  primary 
channel  without  a  signal  installed  to  get  a  measure  of  the 
noise  present  in  the  output.  The  standard  deviation  of  this 
output  of  the  detector  v/nen  this  noise  only  vector  is  us-'d  as 
input  is  used  as  an  approximation  to  the  noise  power. 
Dividing  the  outputs  of  the  detector  in  Figures  11  and  12  by 
this  estimate  of  the  noise  power,  an  estimate  of  the  signal - 
to-noise-ratio  (SNR)  can  be  m.ade.  The  values  of  the  average 
peak  SNR  for  the  three  locations  in  each  run  are  plotted 
versus  signal  strength  in  Figure  13.  The  S!CR  of  the  SFR 
applied  data  runs  ab^^ut  6  dB  higher  than  the  non-SER  applied 


data . 


Figure  9  Demodulator  Output  of  In-Phase  Component; 
(a)  Signal/500  (b)  Signal/1000  (c)  Signal/2000 
(d)  Signal/4000 
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Amplitude  Amplitude  Amplitude 


Figure  10  Demodulator  Output  of  Quadrature  Components; 
(a)  SIGNAL/500  (b)  SIGNAL/1000  (c)  SIGNAL/2000 
(d)  SIGNAL/4000 


51 


Amplitude 


Figure  12  Quadrature  Detector  Outputs:  (a)  Signal/500 
(b)  Signal/1000  (c)  Signal/2000  (d)  Signal/4000 
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Quadrature 


Figure  13  SNR  vs  Relative  Signal  Strength 


D.  DETECTOR  THEORY 

The  process  of  adaptive  noise  cancellation  is  performed 
for  the  purpose  of  increasing  the  signal-to-noise  ratio.  A 
section  on  the  narrowband  detection  schemes  that  have  been 
considered  in  conjunction  with  the  noise  canceler  are 
presented  since  they  are  integral  to  the  analysis  here.  Two 
likelihood  ratio  detectors  are  considered  for  the  detection  of 
narrowband  signals. 
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The  first  detector  is  designed  with  the  a8sun\ption  that 
not  on]y  is  the  signal  known,  but  its  peak  within  a  given 


interval  is  constant  relative  to  the  start  of  that  interval. 
This  last  assumption  is  to  assure  that  the  demodulated  in- 
phase  and  quadrature  signals  that  are  used  in  the  likelihood 
ratio  detector  match  those  in  the  noise  canceled  and 
demiodulated  output  of  the  sensor  system  de.scribed  previously. 
In  the  second  detector,  the  latter  assumption  was  dropped  and 
the  design  followed  that  of  an  optimal  envelope  detector.  In 


both 


ises  It  was  assigned  that  tb.e  residual  noise 


-eld  was 


gaussian.  This  ass’umption,  while  not  strictly  true,  seems  to 
yield  quite  reasonable  results . [Ref .  9] 


The  signal  is  ge.ierated  by  a  calibrated  source  with  knov'm'. 
qualities.  The  modelling  equation  for  the  narrowband  signal 
is  given  by: 


s(t)  =  A  ( t )  cos  ( 2  fp  t  +  6  ( t )  ]  (62) 

For  the  first  detector,  the  phase  term,  0,  is  assurried  constant 
within  any  non-overlapping  interval  0  to  T.  The  interval  of 
0  to  T  is  equivalent  to  the  length  of  the  signal .  The  purpose 
of  this  assumption  is  so  that  a  Matched  Filter  can  be  created 
for  the  in-phase  and  quadrature  components  by  using  the 
demodulated  in-phase  and  quadrature  component  of  the  signal. 
Figure  14  shows  the  demodulated  signals  used  for  the  detector. 

The  first  detector  i?  a  mult -i -channel  Matched  Filter  shown 
in  Figure  15.  The  Matched  Filters,  .h.  and  h,,  are  determaned 


Figure  14  Demodulated  Signal  for  Match  Filtering 


based  on  the  requirement  that  the  output  signal-to-noise  ratio 
is  maximized  at  time  T.  The  corresponding  equations  are 


J  Ki(t-T 


-T)  (t-T)  (t)  _  (T-T) 

-T)  S^(T-t) 


(63) 


where  the  quantities  R  denote  the  various  auto  and  cross 
correlations.  This  detector  is  the  one  used  for  all 
processing  performed  in  this  thesis.  The  input  to  the 
detector  is  the  noise  cancelled  output  of  the  SER  portion  of 
the  processor  after  demodulation . [Ref .  10] 

The  second  detector  is  an  optimal  nonlinear  envelope 
detector.  The  description  of  its  principles  and  design  are 
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Figure 


Noise  Cancelled,  \ 

Demodulated, - 

Quadrature  i 


Noise  Cancelled,  i 
Demodulated,  -- — ht 
In-Phase  I 


Detector 

Output 


15  Linear  Multi “Channel  Match  Filter 


located  in  the  Appendix.  The  detector  was  not  used  by  the 
author  in  v.his  work  relating  to  the  analysis  of  the  SER 
dlgorithrr,,  put  it  is  presented  as  con^pieDentary  work  relating 
to  the  overall  project. 


E.  SIGNAL  LOCATION 

The  signal  strength  analysis  showed  that  the  signal 
divided  by  2000  is  easily  discernable  from  the  background  at 
the  output  of  the  detector.  The  effect  of  varying  the 


loca 

tion 

cf  the  signal  wi 

11  allow 

analysis 

to 

get 

a  better 

feel 

for 

the  expected  SNR 

based  or. 

a  larger 

sar 

r.ple 

size  for 

aver 

aginc 

f  than  previously 

used . 

A  second 

ber. 

,ef  it 

of  this 

analysis  is  to  see  the  effect  of  the  initialization  process  on 
the  SNR.  The  signal  will  be  shifted  by  a  quarter  of  its 
length  for  each  run,  or  2048  points.  The  SNR  is  plotted 
versus  the  location  of  the  center  of  the  signal  in  Figure  16. 

The  data  is  run  with  and  without  the  SER  algorithjn  to  get 


a  good  feel  for  the  enhancement  achieved  with  the  SER 


Figure  16  SI'TR  vs  Location  of  Signal  In  Noise  Vector  (Y) 


algorittjn.  Three  separate  locations  of  the  signal  will  be 
used  for  each  run  to  reduce  calculation  time.  The  use  of 
multiple  signals  causes  no  effect  on  the  SNR  calculation  as 
long  as  the  signals  are  not  allowed  to  overlap. 
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•T  TTT  T'^r  V 


VII.  REStJLTS  A^©  CONCLUSIC^I 


A.  ANALYSIS  OF  RESULTS 

The  use  of  the  SER  algorithm  enhances  the  ability  to 
detect  narrowband  low  power  signals  buried  in  the  background 
noise.  Figure  17  shows  the  spectrum  of  the  noise  sample  from 
sensor  XE  before  and  after  the  use  of  the  SER  algorithm.  The 
frequency  range  is  from  one  to  ten  hertz,  which  is  in  the 
irsriQS  of  conc0irn  fcir  t'no  oHV  s©^n  st 

hertz  (note  the  x-axis  is  a  logarithmic  scale)  is  caused  by 
the  processing  equipment.  This  is  a  man-made  noise  source  and 
will  be  filtered  out  in  the  final  project  design.  The  figure 
illustrates  the  10-15  dB  of  noise  cancellation  possible  by 
using  this  method  of  adaptive  noise  cancellation. 

The  location  of  the  signal  in  comparison  to  the  start  of 
the  adaptive  process  has  been  shown  to  have  very  little  effect 
on  the  signal  enhar.ceme.nt .  The  system  implementation  will 
most  likely  cause  the  algorithm  to  only  require  initialization 
once  a  week  and  operate  constantly  otherwise.  The 
initialization  process  takes  about  109  iterations  or  1.67 
mdnutes  before  settling  out  but  is  heavily  data  dependent. 
Even  while  settling,  the  performance  of  the  algorithm  is 
adequate  to  receive  .mcst  signals.  For  the  purpose  of  this 


project,  initialization  should  be  of  no  concern. 
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Figure  17  illustration  of  SER  Algorithrri  Cancellation 
Capability,  Power  Spectral  Density  of  1-lC  Hz 

B.  COlJCLaSIONS  AND  RECOMMENDATIONf 


It  was  hoped  at  the  beginnir.g  of  the  research  that 
inf orrration  could  be  gained  on  the  characteristics  of  the 
noise  environinent .  The  stationarity  of  the  background  noi.se 
environinent  was  found  to  be  at  least  nine  seconds.  This  is 
be.lievGd  to  be  only  a  characteristic  of  the  specific  data  set 
that  is  used.  The  length  of  stationarity  affects  the  value  of 
the  forgetting  factor,  a,  and  is  implied  by  the  use  of  the  SEP 
algorithm,.  A  better  feel  for  a  good  average  value  to  use  for 
the  length  of  stationarity  will  result  from,  the  use  of  m.ore 
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data  sets.  This  analysis  needs  to  be  conducted  and  could  be 
a  source  of  further  work  in  this  area. 

A  second  area  of  work  that  cculd  develop  from  this 
research  is  a  study  of  the  cancellation  possible  from  a 
vertical  above“Surf ace  sensor  and  a  submerged  horizontal 
sensor.  The  purpose  of  this  would  be  to  differentiate  between 
surface  and  submerged  activity.  This  may  prove  to  be 
impossible  but  it  would  be  interesting  to  see  the  amount  of 
correlation  possible  between  the  two  locations.  This  type  of 
arrangement  coujd  also  lead  to  a  v;ider  distance  between 
sensors  which  could  also  increase  the  knowledge  of  the 
background  enviromr.ent . 

The  purpose  of  the  research  presented  in  this  thesis  is  to 
develop  a  method  to  adaptively  cancel  the  ELF  noise  present  in 
a  submerged  sensor  array.  The  SER  algorithm  appears  to  be  an 
adequate  method  to  accomplish  this  goal  of  noise  cancellation. 
Tiie  cancellation  ratio  in  the  frequency  range  of  interest  is, 


on  average,  12.5  dB ,  This  cancellation  ratio  along  with  the 
use  of  optimal  detectors  will  assist  the  system  in  detecting 
narrowband  signals  buiied  deep  in  the  background  noise. 
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APPENDIX 


The  second  detector  used  in  conjunction  v/ith  this  project 
o£  measuring  the  geoelectric  background  noise  is  an  optimal 
nonlinear  envelope  detector.  The  sy.Tbolism  used  3S  slightly 
different  but  efforts  are  made  to  define  all  relevant  terms. 
Th.e  in-phase  and  quadrature  data  (after  norse  cancellation) 
are  denoted  by  r;  and  r^,  respectively.  Similarly  the  in-phase 
and  quadrature  noises  are  represented  by  n,  and  n, .  Assuming 
a  joint  Gaussian  distribution  we  have  the  following  two 
hypothesis: 

//q  !  (ad  ,  ad  )  =  (  ii‘  ,  n'’  ) 

:  (  ar  ,  ad  )  =  (  n'  n-  )  +  (  a'  ,  s'*  ) 

where  denotes  the  null  hypothe.sis  (noise  only)  and  H- 
denotes  the  alternative  (signal  and  noise  present).  The 
vector  notation  meaiiS 


T  ^  . (rA-,r,^^)] 

where  A-  denotes  the  length  of  the  interval  for  which  the 
decision  is  being  made.  It  is  then  clear  that  the  probability 
density  function  for  the  received  data  is  as  follows: 

p  (r|//.^)  «  exp 

where  denotes  the  covariance  matrix  of  the  noise  vector. 
In  this  particular  case,  it  is  assumed  that  the  in-phase  and 


2  "“nn 


6?. 


the  quadrature  components  have  been  decorrelated  and  whitened 
with  the  same  variance  0^  so 


P  ir 


2e^ 


a 


Introducing  polar  variables  {  rj  ,  )  and  (  ^  )  it  is 

found  thac. 


P  (1*1/4) 


2za 


2N 


20‘ 


K  N  S 

E  “25;  rj,5*cos((t»J 

[A”l  >t"i  *“1 


wheiS  r/  =  +  rj'  and  tan  (  )  =  if  /  rj  .  Integrating 

the  above  over  the  angular' variables  will  give' 


P  (Jfl/^:) 


where  Iq  is  the  modified  bessel  function  of  order  zero.  The 
vector  r  on  the  left  hand  side  now  only  refers  to  the  various 
envelopes,  r  =  [  r. ,  .  .  .,  I'.J  .  The  same  is  true  for  the  signaJ 
envelope.  The  likelihood  ratio  detector  is  found  to  be; 


A  - — ^  ^  =exp| 
P'.i-i^'c) 


The  log-iike'i  inocd  ratio  is  then: 


N 


/ 


logX-=- 


N 


N 


E 


^1 

jf 
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The  first  term  is  a  constant,  dependent  only  on  the  known 
signal,  and  can  be  dropped.  Finally  the  test  statistic  or 
detector  output  is [Ref  10]: 


K 

logX^^  log 


The  implementation  and  analysis  of  this  detector  scheme  is 
left  as  a  possible  follow-on  thesis  topic. 
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